Marine aquaculture has the potential to play an important role in the global food supply as a more sustainable protein option than land-based meat production.1 Gentry et al. mapped the potential for marine aquaculture globally based on multiple constraints, including ship traffic, dissolved oxygen, bottom depth .2
For this assignment, you are tasked with determining which Exclusive
Economic Zones (EEZ) on the West Coast of the US are best suited to
developing marine aquaculture for several species of oysters.
Based on previous research, we know that oysters needs the following
conditions for optimal growth:
We will use average annual sea surface temperature (SST) from the years 2008 to 2012 to characterize the average sea surface temperature within the region. The data we are working with was originally generated from NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1.
To characterize the depth of the ocean we will use the General Bathymetric Chart of the Oceans (GEBCO).3
We will be designating maritime boundaries using Exclusive Economic Zones off of the west coast of US from Marineregions.org.
Below is an outline of the steps you should consider taking to achieve the assignment tasks.
To start, we need to load all necessary data and make sure it has the coordinate reference system.
here packagewc_regions_clean.shp)average_annual_sst_2008.tifaverage_annual_sst_2009.tifaverage_annual_sst_2010.tifaverage_annual_sst_2011.tifaverage_annual_sst_2012.tifdepth.tif)# loading in the libraries
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 0.3.5
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(stars)
## Loading required package: abind
library(tmap)
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
##
## The following object is masked from 'package:dplyr':
##
## select
library(terra)
## terra 1.6.17
##
## Attaching package: 'terra'
##
## The following object is masked from 'package:tidyr':
##
## extract
library(tmap)
library(ggplot2)
library(here)
## here() starts at /Users/colleenmccamy/Documents/MEDS/EDS_223_Spatial_Data/assignments/assignment4-colleenmccamy
library(OpenStreetMap)
library(tmaptools)
# reading in the sea surface temperature data
aa_sst_08 <- rast(here("data", "average_annual_sst_2008.tif"))
aa_sst_09 <- rast(here("data", "average_annual_sst_2009.tif"))
aa_sst_10 <- rast(here("data", "average_annual_sst_2010.tif"))
aa_sst_11 <- rast(here("data", "average_annual_sst_2011.tif"))
aa_sst_12 <- rast(here("data", "average_annual_sst_2012.tif"))
class(aa_sst_08)
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
# combining the sst data
list_sst <- list(aa_sst_08, aa_sst_09, aa_sst_10, aa_sst_11, aa_sst_12)
sst_rast <- rast(list_sst)
# adding in a crs
crs(sst_rast)
## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"unknown\",\n ELLIPSOID[\"WGS84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]"
# plottig to check out the data
tm_shape(sst_rast) +
tm_raster()
# loading in the bathymetry data
depth_rast <- rast(here("data", "depth.tif"))
crs(depth_rast)
## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]"
plot(depth_rast)
#bath_sf <- st_as_sf(depth_rast)
print(paste0("The coordinate refence system for both raster data are both WGS 84."))
## [1] "The coordinate refence system for both raster data are both WGS 84."
crs(sst_rast)
## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"unknown\",\n ELLIPSOID[\"WGS84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]"
crs(depth_rast)
## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]"
Next, we need process the SST and depth data so that they can be combined. In this case the SST and depth data have slightly different resolutions, extents, and positions. We don’t want to change the underlying depth data, so we will need to resample to match the SST data using the nearest neighbor approach.
# finding the mean SST from 2008-2012
sst_mean <- mean(sst_rast)
tm_shape(sst_mean)+
tm_raster() +
tm_layout(main.title = "Mean Sea Surface Temperature from 2008-2012")
# calculating the mean from Kelvin to Celsius
sst_mean_c <- sst_mean - 273.15
# cropping depth data to sst raster extent
depth_rast_crop <- crop(depth_rast, sst_mean_c)
# matching the resultion of both rasters
depth_rast_crop_res <- resample(x = depth_rast_crop,
y = sst_mean_c,
method = "near")
#checking the CRS
crs(sst_mean_c) == crs(depth_rast_crop_res)
## [1] FALSE
# stacking the two rasters of depth and sea surface temperature
list_sst_depth <- list(sst_mean_c, depth_rast_crop_res)
sst_depth_rast <- rast(list_sst_depth)
# plotting the new raster
# plot(sst_depth_rast)
In order to find suitable locations for marine aquaculture, we’ll need to find locations that are suitable in terms of both SST and depth.
1 and unsuitable values to
NAlapp() function
multiplying cell values# creating a reclassification matrix valid sst locations
rcl_sst <- matrix(c(-Inf, 11, NA,
11, 30, 1,
30, Inf, NA),
ncol = 3,
byrow = TRUE)
# using the reclassifying matrix to set non-suitable sst to NA
sst_suit <- classify(sst_mean_c, rcl = rcl_sst)
# creating a reclassification matrix valid depth locations
rcl_depth <- matrix(c(-Inf, -70, NA,
-70, 0, 1,
0, Inf, NA),
ncol = 3,
byrow = TRUE)
# using the reclassifying matrix to set non-suitable depth to NA
depth_suit <- classify(depth_rast_crop_res, rcl = rcl_depth)
# cropping sst and depth data based on the mask
sst_oyster <- crop(sst_mean_c, sst_suit)
depth_oyster <- mask(depth_rast_crop_res, depth_suit)
# reprojecting depth_suit as sst_suit due to error below
depth_suit <- project(depth_suit, crs(sst_suit))
# combining the two rasters
list_oyster <- list(depth_suit, sst_suit)
sst_depth_oyster <- rast(list_oyster)
# overlaying the data
fun_oyster <- function(x, y){return(x * y)}
sst_depth_suitable <- lapp(sst_depth_oyster,
fun = fun_oyster)
#sst_depth_oyster_apply <- sst_depth_oyster$depth * sst_depth_oyster$mean
tm_shape(sst_depth_suitable) +
tm_raster() +
tm_layout("Suitable Locations for Oysters based on Sea Surface Temp. and Depth")
We want to determine the total suitable area within each EEZ in order to rank zones by priority. To do so, we need to find the total area of suitable locations within each EEZ.
#loading in the data and plotting to check it out
wc_eez_vect <- vect(here("data", "wc_regions_clean.shp"))
# plot(wc_eez_vect)
# wc_eez_vect
wc_eez_vect_area <- shapefile(here("data",
"wc_regions_clean.shp"))
wc_eez_vect_area$Area <- area(wc_eez_vect_area)
#print(wc_eez_vect_area$Area)
wc_eez_df <- as_tibble(wc_eez_vect_area) |>
dplyr::select(rgn_key, area_km2)
print(paste("The area of the grid cells are",
wc_eez_df[1, 1, 1], "with an area of",
round(wc_eez_df[1, 2, 1]), "km^2",
wc_eez_df[2, 1, 2], "with an area of",
round(wc_eez_df[2,2,1]), "km^2",
wc_eez_df[3, 1, 2], "with an area of",
round(wc_eez_df[3,2,1]), "km^2",
wc_eez_df[4, 1, 2], "with an area of",
round(wc_eez_df[4,2,1]), "km^2",
"."))
## [1] "The area of the grid cells are OR with an area of 179994 km^2 CA-N with an area of 164379 km^2 CA-C with an area of 202738 km^2 CA-S with an area of 206861 km^2 ."
# converting to a raster
nams <- names(wc_eez_vect)
wc_eez_rast <- lapply(nams, function(x) {
rasterize(wc_eez_vect, sst_depth_suitable,
field = x,
touches = TRUE)})
# merging all objects into one raster
wc_eez_rast <- do.call("c", wc_eez_rast)
#plot(wc_eez_rast)
# using the suitable locations to mask the wc_eez_rast
wc_eez_rast <- project(wc_eez_rast,
sst_depth_suitable)
# finding the area of a raster cell
cell_area <- cellSize(wc_eez_rast, unit = "km")
# plotting cell size
tm_shape(cell_area) +
tm_raster() +
tm_layout(main.title = "Plotting Cell Size")
# masking the suitable areas with the wc_eez_raster
wc_eez_suitable <- mask(wc_eez_rast$rgn,
sst_depth_suitable)
#plot(wc_eez_suitable)
# extracting the areas
area_suitable <- expanse(wc_eez_suitable, unit = "km")
print(paste0("The total suitable area for oysters based on sea surface temperature and depth is ", round(area_suitable), " km^2."))
## [1] "The total suitable area for oysters based on sea surface temperature and depth is 11526 km^2."
# extracting area per suitable region and turning it into a dataframe
area_region <- expanse(wc_eez_suitable, unit = "km", byValue = TRUE)
area_region <- as_tibble(area_region)
# adding a region column in the dataframe
region <- tribble(
~region, ~rgn_key,
"Central California", "CA-C",
"Northern California", "CA-N",
"Oregon", "OR",
"Southern California", "CA-S",
"Washington", "WA")
area_region <- cbind(area_region, region)
# printing the answers
print(paste0("The suitable area for oysters in the ", area_region$region[1], " region is ", round(area_region$area[1]), " km^2."))
## [1] "The suitable area for oysters in the Central California region is 4070 km^2."
print(paste0("The suitable area for oysters in the ", area_region$region[2], " region is ", round(area_region$area[2]), " km^2."))
## [1] "The suitable area for oysters in the Northern California region is 178 km^2."
print(paste0("The suitable area for oysters in the ", area_region$region[3], " region is ", round(area_region$area[3]), " km^2."))
## [1] "The suitable area for oysters in the Oregon region is 1074 km^2."
print(paste0("The suitable area for oysters in the ", area_region$region[4], " region is ", round(area_region$area[4]), " km^2."))
## [1] "The suitable area for oysters in the Southern California region is 3812 km^2."
print(paste0("The suitable area for oysters in the ", area_region$region[5], " region is ", round(area_region$area[5]), " km^2."))
## [1] "The suitable area for oysters in the Washington region is 2393 km^2."
# calculating the percent of each area in the suitable region by joining the two dataframes and then calculating the percent per suitable region
area_per <- left_join(area_region, wc_eez_df, by = "rgn_key")
area_per <- area_per |>
mutate("area_percent" = area/area_km2 * 100)
# printing the results
print(paste0("The percent of suitable area for oysters in the ", area_per$region[1], " region is about ", round(area_per$area_percent[1], 2), "%."))
## [1] "The percent of suitable area for oysters in the Central California region is about 2.01%."
print(paste0("The percent of suitable area for oysters in the ", area_per$region[2], " region is about ", round(area_per$area_percent[2], 2), "%."))
## [1] "The percent of suitable area for oysters in the Northern California region is about 0.11%."
print(paste0("The percent of suitable area for oysters in the ", area_per$region[3], " region is about ", round(area_per$area_percent[3], 2), "%."))
## [1] "The percent of suitable area for oysters in the Oregon region is about 0.6%."
print(paste0("The percent of suitable area for oysters in the ", area_per$region[4], " region is about ", round(area_per$area_percent[4], 2), "%."))
## [1] "The percent of suitable area for oysters in the Southern California region is about 1.84%."
print(paste0("The percent of suitable area for oysters in the ", area_per$region[5], " region is about ", round(area_per$area_percent[5], 2), "%."))
## [1] "The percent of suitable area for oysters in the Washington region is about 3.58%."
Now that we have results, we need to present them!
Create the following maps:
Include:
# combining area and percent area with the vector and setting it as an sf object
map_data <- merge(wc_eez_vect, area_per, by = "rgn_key")
map_data_sf <- st_as_sf(map_data)
suitable_area <- map_data_sf
percent_suitable_area <- map_data_sf
tmap_mode("view")
## tmap mode set to interactive viewing
# setting the basemap
tmap_style("natural")
## tmap style set to "natural"
## other available styles are: "white", "gray", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"
# mapping the data with area and percent of area as the fill
tm_shape(suitable_area) +
tm_fill("area", title = "Total Suitable Area",
palette= (c("#90e0ef", "#48cae4", "#00b4d8",
"#0096c7", "#0077b6", "#023e8a"))) + tm_borders(col = "#000e30")+
tm_layout("Suitable Area for Oysters by EEZ Region" ) +
tm_shape(percent_suitable_area) +
tm_fill("area_percent", title = "Percent of Area Suitable",
palette = (c("#92e6a7", "#4ad66d", "#25a244",
"#208b3a","#1a7431","#155d27"))) + tm_borders(col = "#064208")
print("Use layer toggles to switch between total suitable area and the percent of total suitable area per EEZ region.")
## [1] "Use layer toggles to switch between total suitable area and the percent of total suitable area per EEZ region."
Now that you’ve worked through the solution for one group of species,
let’s update your workflow to work for other species. Please create a
function that would allow you to reproduce your results for other
species. Your function should be able to do the following:
Run your function for a species of your choice! You can find information on species depth and temperature requirements on SeaLifeBase. Remember, we are thinking about the potential for marine aquaculture, so these species should have some reasonable potential for commercial consumption.
https://www.sealifebase.ca/summary/Patiria-miniata.html
print("In order to use this funciton make sure you have read in and processed the data in the first part of this document.")
## [1] "In order to use this funciton make sure you have read in and processed the data in the first part of this document."
fun_suitable_range <- function(temp_min, temp_max, depth_min, depth_max, species){
rcl_sst_sp <- matrix(c(-Inf, temp_min, NA,
temp_min, temp_max, 1,
temp_max, Inf, NA),
ncol = 3,
byrow = TRUE)
sst_suit_sp <- classify(sst_mean_c, rcl_sst_sp)
rcl_depth_sp <- matrix(c(-Inf, depth_max, NA,
depth_max, depth_min, 1,
depth_min, Inf, NA),
ncol = 3,
byrow = TRUE)
depth_suit_sp <- classify(depth_rast_crop_res,
rcl = rcl_depth_sp)
sst_crop_sp <- crop(sst_mean_c, sst_suit)
depth_mask_sp <- mask(depth_rast_crop_res, depth_suit_sp)
depth_suit_sp <- project(depth_suit_sp, crs(sst_suit_sp))
list_suit_sp <- list(depth_suit_sp, sst_suit_sp)
sst_depth_suit_sp <- rast(list_suit_sp)
fun_suitable_sp <- function(x, y){return(x * y)}
sst_depth_suitable_sp <- lapp(sst_depth_suit_sp,
fun = fun_suitable_sp)
wc_eez_suitable_sp <- mask(wc_eez_rast$rgn,
sst_depth_suitable_sp)
area_suitable_sp <- expanse(wc_eez_suitable_sp, unit = "km")
print(paste0("The total suitable area for ", species, " based on sea surface temperature and depth is ",
round(area_suitable_sp), " km^2."))
area_region_sp <- expanse(wc_eez_suitable_sp, unit = "km", byValue = TRUE)
area_region_sp <- as_tibble(area_region_sp)
area_region_sp <- cbind(area_region_sp, region)
print(paste0("The suitable area for ", species, " in the ", area_region_sp$region[1], " region is ",
round(area_region_sp$area[1]), " km^2."))
print(paste0("The suitable area for ", species, " in the ", area_region_sp$region[2], " region is ", round(area_region_sp$area[2]), " km^2."))
print(paste0("The suitable area for ", species, " in the ", area_region_sp$region[3], " region is ",
round(area_region_sp$area[3]), " km^2."))
print(paste0("The suitable area for ",species, " in the ", area_region_sp$region[4], " region is ",
round(area_region_sp$area[4]), " km^2."))
print(paste0("The suitable area for ", species, " in the ", area_region$region[5], " region is ",
round(area_region$area[5]), " km^2."))
area_per_sp <- left_join(area_region_sp, wc_eez_df, by = "rgn_key")
area_per_sp <- area_per_sp |>
mutate("area_percent" = area/area_km2 * 100)
print(paste0("The suitable area for ", species, " in the ", area_region_sp$region[1], " region is ",
round(area_per_sp$area_percent[1], 2), "%"))
print(paste0("The suitable area for ", species, " in the ", area_region_sp$region[2], " region is ",
round(area_per_sp$area_percent[2], 2), "%"))
print(paste0("The suitable area for ", species, " in the ", area_region_sp$region[3], " region is ",
round(area_per_sp$area_percent[3], 2), "%"))
print(paste0("The suitable area for ", species, " in the ", area_region_sp$region[4], " region is ",
round(area_per_sp$area_percent[4], 2), "%"))
print(paste0("The suitable area for ", species, " in the ", area_region_sp$region[5], " region is ",
round(area_per_sp$area_percent[5], 2), "%"))
map_data_sp <- merge(wc_eez_vect, area_per_sp, by = "rgn_key")
map_data_sp_sf <- st_as_sf(map_data_sp)
suitable_area_sp <- map_data_sp_sf
percent_suitable_area_sp <- map_data_sp_sf
tmap_mode("view")
tmap_style("natural")
print(paste0("Suitable Area for ", species, " by EEZ Region"))
tm_shape(suitable_area_sp) +
tm_fill("area", title = "Total Suitable Area",
palette = (c("#90e0ef", "#48cae4", "#00b4d8",
"#0096c7", "#0077b6", "#023e8a"))) + tm_borders(col = "#000e30")+
tm_shape(percent_suitable_area_sp) +
tm_fill("area_percent", title = "Percent of Area Suitable",
palette = (c("#92e6a7", "#4ad66d", "#25a244",
"#208b3a", "#1a7431", "#155d27")))+ tm_borders(col = "#064208")
#print("Use layer toggles to switch between total suitable area and the percent of total suitable area per EEZ region.")
}
# The bat star has a sea surface temperature range of approximately 2-20 degrees (C) and can live in depths from 0-302 feet below sea level. They are an interesting species to look at as they are affected by sea star wasting disease which is correlated to an increase in ocean temperatures.
fun_suitable_range(temp_min = 2, temp_max = 20,
depth_min = 0, depth_max = -302,
species = "Bat Star")
## [1] "The total suitable area for Bat Star based on sea surface temperature and depth is 67929 km^2."
## [1] "The suitable area for Bat Star in the Central California region is 10322 km^2."
## [1] "The suitable area for Bat Star in the Northern California region is 9029 km^2."
## [1] "The suitable area for Bat Star in the Oregon region is 18146 km^2."
## [1] "The suitable area for Bat Star in the Southern California region is 12292 km^2."
## [1] "The suitable area for Bat Star in the Washington region is 2393 km^2."
## [1] "The suitable area for Bat Star in the Central California region is 5.09%"
## [1] "The suitable area for Bat Star in the Northern California region is 5.49%"
## [1] "The suitable area for Bat Star in the Oregon region is 10.08%"
## [1] "The suitable area for Bat Star in the Southern California region is 5.94%"
## [1] "The suitable area for Bat Star in the Washington region is 27.12%"
## tmap mode set to interactive viewing
## tmap style set to "natural"
## other available styles are: "white", "gray", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"
## [1] "Suitable Area for Bat Star by EEZ Region"
Hall, S. J., Delaporte, A., Phillips, M. J., Beveridge, M. & O’Keefe, M. Blue Frontiers: Managing the Environmental Costs of Aquaculture (The WorldFish Center, Penang, Malaysia, 2011).↩︎
Gentry, R. R., Froehlich, H. E., Grimm, D., Kareiva, P., Parke, M., Rust, M., Gaines, S. D., & Halpern, B. S. Mapping the global potential for marine aquaculture. Nature Ecology & Evolution, 1, 1317-1324 (2017).↩︎
GEBCO Compilation Group (2022) GEBCO_2022 Grid (doi:10.5285/e0f0bb80-ab44-2739-e053-6c86abc0289c).↩︎